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ABSTRACT 

We present a detailed study of the collapse of molecular cloud cores using high resolution 
3D adaptive mesh refinement (AMR) numerical simulations. In this first in a series of inves- 
tigations our initial conditions consists of spherical molecular core obeying the hydrostatic 
Bonnor-Ebert-Profile with varying degrees of initial rotation. Our simulations cover both the 
formation of massive disks in which massive stars form as well as low mass disks. We use a 
customized version of the FLASH code whose AMR technique allows us to follow the for- 
mation of a protstellar disk and protostellar core(s) through more than ten orders in density 
incre ase w hile continuously resolving the local Jeans length (i.e. obeying the Truelove crite- 
rion, l^raeloveetalj 119231)) ■ Our numerical simulations also incorporate the energy loss due 
to molecular line emission in order to obtain a more realistic picture of protostellar core and 
disk formation. Out initial states model system of mass 168 Mq and 2.1 Mq that will form 
high and low mass stars, respectively. We follow many features such as the development com- 
plex shock structures, and the possible fragmentation of the disk. We find that slowly rotating 
cores (fl tff — 0.1) produce disks in which a strong bar develops but which does not fragment. 
Faster initial rotation rates (ft tff = 0.2) result in the formation of a ring which may fragment 
into two star-forming cores. The size of the rings found in our simulated disks agree with the 
observations of similar systems. 

Key words: accretion, accretion discs, hydrodynamics, ISM: clouds, evolution, methods: 
numerical 



1 INTRODUCTION 

The formation of protostars and their surrounding disks from 
collapsing, dense, molecular cores within molecular clouds 
is still an actively debated top i c. Early pion eerin g works 
by Bode nheimer & SweigartI il968h . iLarsonl |l969'), and lPenstonI 
U969..) addressed this problem using numerical sim ulations for var- 
ious initial conditions. iBodenheimer & Sweigarl studied collapses 
of isothermal spherical symmetric gas clouds which are initially 
out of equilibrium. They pointed out that the solutions depend only 
very weakly on the choice of the initi al density profile (homoge- 
neous or Chandrasekhar-tvpe). iLarsorl il969^ and .Penston U969.) 
included the effect of energy dissipation by radiation in their ID 
simulations and found that cooling during the initial collapse phase 
is efficient enough to keep the gas at its initial temperat ure. Al- 
though both studies start with different initial cond itions jLarsonI 
used a homogeneous density distribution whereas IPenstonI used 
a Bonnor-Ebert-Profile) they obtained similar results during the 
isothermal collapse. Starting with a highly unstable, initial singular 



isothermal sphere IShd il977h obtained a different evolution sce- 
nario: a self-similar inside-out collaps e and c l aimed that the for- 
mer solutions are physically artificial. iHunteJ ^9^) in his study 
of unstable isothermal spheres pointed out that both investigations 
address the same proble m but refer to two different stages of the 
isothermal collapse. The iLarsonllPenstonl studies investigate the 
collapse of molecular clouds prior to the formation of protostel- 
lar core whereas Shu's study apply to the accretion phase of an 
non-rotating cloud with a central protostar 

In reality, molecular cores have some inital angular momen- 
tum so that most of the collapsing material ends up in a rotating ac- 
cretion disk (due to their non-vanishing initial angular momentum). 
This significantly alters the dynamics of the protostellar core's evo- 
lution. More recently, this issue was addressed with the help of 2D 
and 3D numerical s imulations by several groups (for a review see 
IBodenheimer et all I2OOO1) . For instance iBurkert & BodenheimeJ 
^^93ir studied the fragmentation of a rotating sphere of uniform 
density with a m — 2 perturbation on a fixed nested grid. Their 
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simulations showed thiat a bar formed whicii tiien fragmented into 
two low mass binari es and several lighter secondary fragments. 
lTrueloveetai] h998'^ used a 3D AMR technique that resolves the 
local Jeans length throughout the simulation to reinvestigate the 
fragmentation problem. These authors found that a bar like struc- 
ture collapses instead to a single filament without further fragmen- 
tation. These differences can be shown to depend upon the ability 
to resolve the local Jeans length during the collapse (e.g. review 
[Sodenheimer et al. 2000). An equivalent requ irement for SPH sim- 
ulations was found by B ate & Burkerd ^993) who showed that the 
minimal resolvable mass must be smaller than the Jeans mass to 
avoid numerical fragmentation. But these authors also noted that 
the inclusion of artificial viscosity and heating, which increases the 
Jeans length and slows the collapse, lead to physical fragmenta- 
tion of bars which is predicted in the work of llnutsuka & Mivamal 
il99 2h. Clouds with an initial Gaussian density profile were studied 
by [Boss 1 1993) who concluded that the most sensitive parameter 
for clouds to fragment is a, the ratio of the thermal to gravitational 
ener gy. 

iMatsumoto & Hanawal i2003h did an extensive parameter 
study of rotating Bonnor-Ebert-Spheres with different rotation pro- 
files and rotation speeds. Their investigations are based on a three 
dimensional nested grid technique which resolves the local Jeans 
length and the use of two different equation of states depending on 
the central density (an isothermal EOS in the low density regime 
and a barotropic EOS in the high density regime). They showed 
that the final structure of the collapsed cloud (bar, ring, or binary 
system) depend sensitively on the parameter Uf Q., where tn and Q, 
are the initial free fall time and the angular velocity, respectively. 

Another important physical process that alters the dynamics 
of the core formation and disk evolution is the cooling process dur- 
ing the collapsing phase. An often used simplification to account 
for the decreasing efficiency of the cooling mechanism is the previ- 
ously mentioned sudden switch of the EOS, depending on the local 
density. We do not follow this approach but rather include the loss 
of energy due to molecular line emission by coUisional excitations. 
In this work we show that the effective equation is a complex func- 
tion of time, density, and space which influences the structure of 
the protoplanetary disk and its possible fragmentation. 

The study of Bonnor-Ebert-Spheres as initial states for star 
formation are of particular interest because they resemble marginal 
stable molecular cores in pressure equilibrium and have flat 
topped rather than singular density profiles. Observed evidence for 



Bonnor-Ebert-type cores are found in a variety of molecular clouds 


(e.g. 


Racca et al. 


2002: Harvev et al. 2001: Alves et al. 


2001,: 


) and 



in num erical simulations of star cluster formation in molecular 
clouds frillev & Pudritzl2004 . The oberservation of massive disks 
undergoing high mass star formation iChini et alj E004) suggests 
that a similar scenario pertains to massive star formation (e.g. 
lYorke & Soiinhalter 2002) . 

In this paper we investigate low and high mass star forma- 
tion from rotating molecular cores through AMR simulations. This 
paper is organized as follows: in Sec.|2|we explain our numerical 
scheme including the cooling procedure, and summarize the prop- 
erties of the static Bonnor-Ebert-Sphere. Our results of the collapse 
in the isothermal and non-isothermal regime are given in Sec.|3| 
This sections includes also a review of the pure isothermal collapse 
of a overcritial Bonnor-Ebert-Sphere. The formation of disks in low 
and high mass rotating cores are discussed in Sec.|4]and in Sec.|5| 
we present our results on the formation of rings and bars in the 
disk and their possible fragmentation. Finally, we summarize and 
discuss our findings in Sec.|6| 



2 NUMERICAL SCHEME AND INITIAL CONDITIONS 

For o ur studies of collapsing gas clouds we used the FLASH 
code jFrvxell et all [2000) which solves the coupled gravito- 
hydrodynamic equations on an adaptive mesh. FLASH is based 
on a block structured adaptive mesh refinement (AMR) technique 
which is implemented in the PARAMESH library iOlson et aU 
iI999i) . This AMR technique allows us to follow the core formation 
over more the ten orders of magni tude in density incre ase with- 
out violating the Truelove criterion iTruelove et alll997l) . For this 
purpose we implemented a new refinement criterion to the FLASH 
code which assures that the local Jeans length 
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is at least resolved by A'^ > 4 grid points, where, c, G, and p are the 
isothermal sound speed, the gravitational constant, and the local 
mass density, respectively. The runs we present in this work are 
performed with an even smaller Jeans number J = Ax/Aj of 1/8 
or 1/12 (Aa: is the grid spacing in one dimension at the point x). 

2.1 Properties of Bonnor-Ebert-Spheres 

We choose as an initial setup a non magnetized spherical gas cloud 
in a ma rginal stable hydrostatic equilibrium, i.e. a Bonnor-Ebert- 
Sphere jEberll9.5.5l : lBonnoJl95' 61 With the definitions 
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the radial density profile is given by the solution of the Lane-Emden 
equation (cf. Chandrasekhar 1967.) : 
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where po is the initial central density, c is the thermal sound speed, 
ro is the characteristic radius of the gas cloud, and a prime (') de- 
notes the differentiation with respect to f . The solution is charac- 
terized by a density profile with a flat core of density po and an 
envelope which falls off with radius roughly as r"^ until truncated 
by the pressure of the surrounding medium. Note that the function 
0(C) is - up to an integration constant - equal to the gravitational 
potential ip (cf. Eq. Jilt ). 

We use a fourth order Runge-Kutta numerical scheme to solve 
Eq. J4j for the potential (f){^), and the acceleration (!>'{£,) which 
give the density profile p(f) and the mass parameter q(^). Fig.Q 
summarizes the radial dependences of the dimensionless parame- 
ters which determine the Bonnor-Ebert-Sphere with the cut-off ra- 
dius of ~ 6.5. 

The total mass within a Bonnor-Ebert sphere at the radius ^ is 
given by 
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where q is the integrated dimensionless mass factor (see also 
iMcLauehlin & Pudritzligg^) . 
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Figure 1. Shows the dimensionless parameters of a Bonnor-Ebert-Sphere: 
the density profile (solid line) multiplied by a factor of 10, the potential <j> 
(dashed line), the gravitational acceleration 0' (thin dashed line), the mass 
parameter q (dotted line), and the ratio of the thermal energy to gravitational 
energy a (dash-dotted line) multiplied by a factor of 10 for a sphere with 
= 6.5. The vertical line indicates this cut-off radius. 
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Such a sphere is u nstable to coll apse if its dimensionless radius ^ is 
larger than 6.451 l'Bonnor'l95^^. Given the external pressure Pext 
and the sound speed of such a spherical cloud the physical radius 
and total mass are, respectively: 



P = 49 — P"^/^ 



M 
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The ratio of thermal to gravitational energy, a, within a 
Bonnor-Ebert-Sphere, which varies only slightly with radius, is 
completely fixed by its truncation radius ^c- This can seen by the 
following argument. The gravitational potential V" inside the sphere 
is determined by 



r2 



(10) 



d r 

where the mass M is given by Eq. 0. Using the matching con- 
ditions at the edge of the sphere with the potential in the ambient 
media, i.e. ip oc r^^ for ^ > ^c, the solution of Eq. <10t is 

iP = ~c [<^c + gc - 0] for C < , (11) 

where <^c = <t>{£,c), (t>'c = ^'(^c), and qc = g(Cc)- Therefore, 
a — 3/2 c^/|i/'i is given by 

a = . . (12) 

With ~ 6.5 (</>c ~ 2.66, qc = 15.85) a is in the range of: 

0.081 < « < 0.094 . (13) 

The rotational to gravitational energy /? for a rigidly rotating 
Boimor-Ebert-Sphere is given by 
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^ Note, as 5 = 2nr/Xj this instability ciiterion is almost equal to Jeans' 
instability criteiion, i.e. r > Xj, and the amount of Jeans masses within the 
sphere is S/tt"^ q ~ 1.53. 
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0.2 


1.15 X 10" 


-2 


A3 


3.63 X 10" 


8.265 X 10- 


15 
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2.12 X 10^2 
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Table 1. The initial parameters for the simulations we present in this paper 
Where f2, tff, and /3max are the initial angular velocity, the initial free fall 
time, and the maximum of ratio of the rotational and gravitational energy, 
respectively. 



where and are the toroidal velocity and angular velocity, re- 
spectively. For instance, a sphere rotating with tff Q, = 0.1 has val- 
ues of /3 that span the range 



0<P< 0.0029 



(15) 



and scales roughly proportional to ^ , while (3 takes its max- 
imum at the edge of the sphere and is given by /3max = 

0.11 {ipx' (tfi^f- 



2.2 Initial conditions: low and high mass cores 

In this paper we present the results of four simulations aiming to 
study both low and high mass pre-stellar cores. We start each simu- 
lation with an overcritical Bonnor-Ebert-Sphere which rotates ini- 
tially with a constant angular velocity Q. Three of these initial mod- 
els have the same initial core density po but different initial angular 
velocities (run Al - A3) and a mass of 168 Mq characteristic of 
massive star formation <Chini et alj2'004l) . and one low mass model 
starting with a core density accord ing to the observed Bok globule 
' 3) with a mass of 2.1 Mq and a 



Barnard 68 by Alves et alJ {2 
medium angular velocity (run B68). The initial parameters of these 
runs are summarized in tableQ 

One particularly interesting object for low mass star forma- 
tio n is the Bok globu le Barnard 68 which was extensively studied 
bv lAlves et alj ^2001). The density profile of this cloud exhibits a 
close to perfect Bonnor-Ebert profile with ~ 6.9 ± 0.2 corre- 
sponding to a physical radius R = 1.25 x 10* AU, a core density 
of Po = 1.0 X 10-"gcm-^, a total mass of Mb68 = 2.1 Mq, 
and a temperature of T = 16 K (we set up this run with an exter- 
nal pressure of Pcxt/fca = 2.4 x lO'"" Kcm~^). This corresponds 
to model marked by a diamond shown in Fig.|2| One can see that 
this observed system is nearly immediately scalable from the initial 
state that we use defined by the triangle. To simulate the collapse of 
a realistic molecular cloud we initializ e run B68 with the physical 
parameters given above. According to lLada et alJ i20()3l) Barnard 
68 might rotate slightly with a ratio of the rotational to gravitational 
energy of a few percent. We initialize the sphere with tfffl = 0.2 
corresponding to /3„,ax ~ 0.01. 

Most of our simulations are designed to study the formation 
of massive accretion disks in which massive (O, B) stars may form. 
The collapse of massiv e (up to 120 M^), initially singular cores 
has been simulated by lYorke & SonnhalteJ i2002h who included 
a careful treatment of radiative transfer elfects. Recent observa- 
tions by Chini et al. (2004) have discovered a massive (20 Mq) 
star that is being formed within a very massive, large accretion 
disk (at least 100 Mq of gas) in the young star forming regin of 
M 17. We therefore chose the initial setup for the study of mas- 
sive star and disk formation corresponding to runs Al - A3 as 
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Figure 2. Scaling of the physical parameters of critical Bonnor-Ebert- 
Spheres. The solid line shows the parameters for sphere with T = 20 K, 
and the upper and lower dashed lines show the relations for a sphere with 
T = 10 K and T = 30 K, respectively (Note that the axis labeling for Pcxt 
is only valid for T = 20 K and has to be rescaled for different tempera- 
tures). The triangle (A) refers to our model parameters for Al - A3, th e 
square (□) shows the Coalsack globule 2 observed by Racca et a^ i2002ft . 
and the diamond (o) marks the properties of Barnard 68 (our model B68) 
as given in .Alves et al. i2()0U . 



follows: the cut-off radius is = 6.5 corresponding to a phys- 
ical radius of i? = 1.62 pc {R = 3.34 x 10^ AU), a unper- 
turbed core density po = 3.35 x lO"^'^ gcm~^, and sound speed 
c = 0.408 km s'^ (T = 20 K). This parameters with a 10% over- 
density result in a total mass of the gas cloud of M — 168 Mq, 
and an external pressure Pext/fcfl = 3.2 X 10^ K cm . The inital 
free fall time tff = -y/ 37r/32 G po of our Bonnor-Ebert sphere is 
1.1 X 10^ years'^. The edge of the sphere is defined by a density de- 
crease of the ambient medium by a factor of 100 whereas the pres- 
sure is continuous at the edge of the sphere. Therefore, the ambient 
low density gas is a hundred times warmer (Tl,mb ~ 2000 K) than 
the cold gas cloud and the sound speed there is Camb = 4.08 kms~^. 
To make sure that the evolution of the sphere is not influenced by 
the simulation box boundaries we choose the simulation volume to 
be (lOOro)^ ~ (15 i?)^ ~ (25 pc)^ 

Fig.|2|shows the scaling relations in the isothermal regime for 
the physical mass, radius, core density, and external pressure for 
critical Bonnor-Ebert-Spheres. All spheres at a given temperature 
with the properties indicated by the shown lines are equivalent. 

Radiative cooling by molecular excitation lines is very ef- 
ficient in the low density regime. In fact, molecular clouds will 
maintain the same temperature during their collapse until the den- 
sity rises above n ~ lO^cm"^. We wish to study this transi- 
tion to less efficient cooling very carefully and therefore start the 
collapse of Bonnor-Ebert-Spheres in this low-density, isothermal 
phase. Therefore, we choose the initial density for the runs Al - 
A3 (no ~ 10'^ cm"'') in order t o sample the full ran^ e of densi- 
ties given by the cooling data of iNeufeld & Ka ufman and 
iNeufeld et al. ( 1995 ) . Their cooling data span a number density of 
10^ < n{H2) < 10^" cm"''. We emphasize however that these 
initial physical conditions can be rescaled only as long as the gas 
stays at the same temperature, i.e. contracts isothermally (cf. dis- 



^ Note that the more accurate dynamical time scale for a Bonnor-Ebert- 
Sphere is to = 1 /V^tt G po 0.52 % 



cussion in section |3)- Therefore, we can apply our results to ob- 
served molecular clouds with different initial physical parameters. 

All runs are set up with a 10% m = 2 density perturbation on 
top of a slightly enhanced Bonnor-Ebert-Profile, i.e., 

P = Pbe(1.1 + 0.1 cos(2(^)) , (16) 

where pbe obeys Eq. J4j and ip is the azimuthal angle (the 10% 
overall density enhancement guarantees the collapse of a critical 
BE-sphere since it overwhelms the relatively small amount of rota- 
tion). All spheres rotate initially with a rigid body rotation profile 
with the amount of rotation given in tableQ 

2.3 Radiative cooling 

Earlier (for a review seelsydenheimer et all2OO0l) and recent (e.g. 
iMatsumoto & Hanawall20(}3[) simulations of protostellar disk for- 
mation and fragmentation typically use two types of equations of 
state: (i) an isothermal equation of state throughout, which assumes 
that radiative cooling is always efficient enough to keep the gas 
at the same initial temperature, or (ii) an isothermal EOS that is 
switched to an adiabatic one at sufficiently high density. We did 
some simulations of the latter approach and found problems not 
only due to numerical artifacts but also due to physical difficulties 
because the sudden change of the equation of state violates energy 
conservation (an isothermal EOS corresponds to an infinite heat 
bath whereas an adiabatic EOS assumes a finite thermal energy). 

To achieve a more realistic, and more physically correct 
picture of the formation of protostellar disks, we account in 
our simulations for the effect of cooling by collisional exci- 
tations of gas molecule s. We use cooling functions provided 
by Neufeld & Kaufm^ <1993h : Neufeld et al. 1 1995). These au- 
thors computed the radiative cooling rates as a function of the gas 
temperature, density and optical depth for the temperature and den- 
sity range of 10 <T < 2500 K and 10^ < n{H2) < lO" cm"^ 
Using steady state molecular abundances for the most important 
coolants in molecular clouds (H2, H2O, CO, O2, HCl, C, and O) 
they provide a total cooling rate which can be used to calculate the 
energy loss due to radiative cooling. We use their optical depth pa- 
rameter for an singular isothermal sphere (which is appropriate for 
a sphere with a 1/r^ density profile) 

iVsis = 5.1 X lO^M 5- cm"^ per kms"^ (17) 

V cm / 

to read off the total cooling rate fro m the provided cooling data 
base. The cooling rates provided by INeufeld & Kaufm'ai] il993h 
and Neufeld et al. (1995) are available for densities n{H2) < 
10^" cm^"*. Because the cooling power per H2 molecule in the 
high temperature (T > 100 K) and high density (n > 10^° cm~^) 
regime is nearly independent of the density, we extrapolate their 
data in the high density range assuming that the cooling power is 
only a function of temperature. As the maximum temperatures in 
our simulations are well below the dissociation temperature of H2 
(Tdiss ~ 2000 K) we do not account for this cooling process in the 
present simulations. We also did not consider cooling processes due 
to gas-grain interactions which we will implement in future studies. 
Howev er, our results are similar in character to those of Yorke et all 
i 19951) who did include dust cooling. 

The treatment of the thermal cooling of the gas allows one to 
calculate the internal energy density, e, at any position and time. 
The pressure, P, at each time step is then calculated by: 

P = e(l-7) , (18) 
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Figure 3. Radial profiles of density (first panel) and the radial velocity (sec- 
ond panel) at different times in the isothermal regime. The sphere starts 
to collapse from outside-in leaving al/r^ density profile in the envelope 
while the flat core region shrinks. The maximal infalling velocity appears 
close to the edge of the core and is moving inward with increasing density. 
The shown lines con'espond to (from bottom to top of the first panel, and top 
to bottom for the second panel) t/to = 3.7, 5.6, 7.13, 7.33, 7.36, 7.37, 
where to = l/\/47r G po ~ 6.0 X 10^ years and the data are taken from 
run Al. 



where 7 is the adiabatic index. Throughout this study we use a con- 
stant 7 of 5/3. Including cooling in our simulations the thermal 
energy is not only altered by compressional heating during the col- 
lapsing phase but also by the loss of energy due to radiative emis- 
sion. We found it is useful to compare our approach which uses 
full molecular cooling, with the earlier adiabatic model studies by 
calculating an effective EOS, i.e. 7cff = dP/dp, at each time step. 



3 CORE FORMATION AND EVOLUTION 

3.1 Pure isothermal collapse of Bonnor-Ebert-Spheres 

Early theoretical work on star formation (e.g. lHavashi & Nakand 
1 1965 ) already pointed out that cooling during the initial stage 
of collapsing molecular cloud cores is efficient enough to ensure 
constant temperatures of the gas over 4-5 orders of magnitude 
in density increase. Hence, the isothermal collapse proceeds as 
a free falling contraction of the gas. Our simulations - includ- 
ing molecular cooling - resemble this early isothermal phase until 
n ~ 10^'^ cm~^. Fig.|5|and Fig.Qfrom our simulation Al clearly 
show the evolution in the isothermal phase at low densities. Be- 



cause the solutions of collapsing isothermal Bonnor-Ebert-Spheres 
are rather general we devote the following paragraph to review the 
collapse of a pure isothermal cloud. 

A detailed ID numerical investiga tion of a collapsing i sother - 
mal Bonnor-Ebert-Sphere was done bv lFoster & ChevalieJ il99l) . 
The results of our 3D simulations in the isothermal regime are very 
similar to the results found by Foster & Chevalier ( 1993). Fig.|3| 
shows the time evolution of the density profile and the radial ve- 
locity profile which have the same trend as the ID results: after a 
sound wave propagates through the sphere at a sound crossing time 
of fsc = i I \/4:Tr G po (~ 3.8 X lO'' years in the case of runs An 
and ~ 2.4 x 10^ years in the case of B68) the cloud starts to col- 
lapse from outside-in. Initially the gravitational acceleration of the 
Bonnor-Ebert-Sphere 



GM 



(19) 



is largest at ^ = 3.0, (0'(3.O) = 0.517), and decreases only 
slightly with radius ((/)'(6.5) = 0.375). Soon after the low den- 
sity material from the outer part of the cloud is accelerated towards 
the center the density profile of the sphere starts to change: a high 
density flat core with a envelope builds up. The l/r*^ law 
in the envelope reflects the fact that the sphere maintains a state 
that is close to hydrostatic equilibrium in the envelope. Now the 
gravitational acceleration is highest at the edge of the flat core and 
decreases with 1 jr. Therefore, material at the core boundary is con- 
stantly accelerated whereas material in the envelope stays nearly at 
a constant infalling velocity because the acceleration becomes very 
low in this region. As the core density increases the core shrinks 
continuously and the envelope profile becomes more domi- 
nant. The maximum of the radial velocity appears around the edge 
of the flat core region and is moves inward with time. Although, 
the infalling material becomes supersonic, no shock occurs in this 
isothermal collapsing phase since the central gravitational force ac- 
celerates the gas in front of the supersonic material fast enough to 
prevent a velocity discontinuity. The radial velocity stays close to 
zero at the outer edge of the sphere and at the center of the cloud. 

We emphasize that the collapse of an over-critical Bo nnor- 
Ebert- Sphere differs from the Expansion- Wave solution of IShd 
since the latter is an inside-out collapse of an initially sin- 
gular isothermal sphere (SIS). The collapse of a SIS proceeds in 
a self-similar fashion leading to a velocity profile that incr eases 
with decreasing radial distance to the center (cf. Fig. 2 of ISh j 
1977) whereas the velocity of a collapsing Bonnor-Ebert-Sphere 
possesses always a finite maximum and goes to zero at the center 
of the sphere. 

To get a quantitative estimate of the isothermal collapse de- 
scribed above we approximate the density profile shortly after the 
initial collapse by: 



p{t) 



pcoi'e 

(r/r 

COl 



r > rcoie 



(20) 



where pcore is the core density and rcore is the radial distance to 
the core edge. This approximation is valid as long as the collapse 
is isothermal. The core becomes non-isothermal after the density 
increased above the critical density of n ~ lO'^'^ cm^'' (see sec- 
tionlT2l. 

Using conservation of the mass of the sphere M we can es- 
timate the core radius for a given core density (with the criti- 
cal density of pcore = 10~'^®gcm~^ and a sound velocity of 
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Figure 4. Comparison of the local free-fall time tff with the cooling time 
tcool in different temperature regimes. For a 20 K cloud (thick line) the 
cooling time becomes longer than the free-fall when its number density ex- 
ceeds n ft! 10^'^ cm^'^. T he cooling data are taken from Neuteld et aL 
1 19951) and Neufeld & Kaufm^ Il993h . The temperature range is T = 
10 K - 10^"'' K (from top to bottom) in steps of 0.1 dex. 



c = 0.41 kms"^): 

« J-^^ = ^=^{i.4>'.Y'' , (21) 

V 47rpco,e^ V47r G Pcore 

~ 4.5 X 10^ AU 

where R is the radius of the sphere and we assumed rcmc ^ -R- 
Using this approximations the total core mass is given by 

^ 47r 3 1 rcorc 

Mco,-e « — r,ore Pcore = 3 

= \ , "\ i^^^'^f' (22) 

~ 0.1 Mq . 

3.2 Collapse with molecular cooling: initial isothermal phase 

Theoretical models of collapsing protostars p redict that the forma- 
tion of a protostar occurs in several stages (e.g. lLarsonll969ll2003h . 
Already from Fig. |4] which compares the cooling time scale tcooi 
and the local free fall time t[{, one can infer that a gravitationally 
unstable low density cloud will undergo a isothermal collapse until 
the cooling by molecular line emission becomes inefficient when 
tcooi > tff- For instance, a 20 K cloud cannot be efficiently cooled if 
its core density exceeds more than 10'^ '^ cm"''. After the core den- 
sity reaches this critical value, the high density region of the gas 
cloud contracts almost adiabatically on a much longer time scale. 
Another prediction one can infer from Fig. |4|is the evolution tra- 
jectory of the core in the temperature-density-plane which follows 
the track where t[[ ~ icooi in the inefficient cooling regime. Our 
simulations are in good agreements with this prediction. 

Initially, due to effective cooling, the unstable sphere main- 
tains its initial temperature while collapsing on its local free fall 
time. Fig. |5| shows the evolution trajectory of the cloud's core in 
the temperature-density plane from our simulation Al. The inital 
contraction stage clearly indicates the regime where cooling is fast 
enough to stay at the clouds initial temperature while core den- 
sity increases. A quantification of this stage is shown in Figs. |6| 
and where we compute the effective equation of state, i.e. 
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Figure 5. Evolution trajectory of the core forming region from our simula- 
tion in the temperature-density plane. Initially, the core undergoes a rapid 
isothermal contraction until the cooling by molecular lines becomes ineffi- 
cient, i.e. tjcoi > iff, at a core density ricore ^ cm~^ (from run Al). 
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Figure 6. Shows the thermal pressure, P = P(ncorc), in the core forming 
region. The inital isothermal collapse (indicated by the solid line, p oc n) 
is followed by a non-isothermal contraction with a non-linear pressure re- 
sponse (from run Al). 

7cff = dlogp/dlogp. We calculate the trajectory plots by aver- 
aging the variables in the region where nmax/S < n < nmax, where 
nmax(t) is the maximal density in our simulation at the time t where 
rimaxit) increases continuously with time. 

From Fig.|4]and from our simulations (Fig.|5|and Fig.0 we 
get a core density (i.e. of material interior to radius n-ore) at the 
end of the isothermal phase of nco,c ~ 10^'^ cm~^ and pcmc ~ 
lO^^*" gcm^'^. Using Eq. <21> and Eq. J22> gives a core radius of 
Tcorc ~ 450AU and a core mass of A/core ~ 0.1 Mq. Note that 
these values are independent of the initial mass M and radius R 
as long as the inital sphere obeys the Bonnor-Ebert profile given 
by Eq. J4}. Other investi gations o f spherical collapses of molecular 
clouds (for a review see lLarsoni2003 ) also predict a first 'hydro- 
static core' when the density reaches ~ 10"^" gcm^^ where the 
core mass is almost independent of the clouds initial mass or ra- 
dius. The predicted 'hydrostatic core' mass is about 0.01 Mq. Our 
findings give a larger core mass at the end of the isothermal col- 
lapsing phase as cooling becomes inefficient at lower densities. 

Comparing the value of the core mass at the end of the isother- 
mal phase to the local leans mass, 

pcore ~ 
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Figure 7. The effective equation of state, 7eff = d log p/d log p, in the core 
forming region as a function of the core density (from run Al). 



0.56 Mq, (which is independent of the sphere parameters) gives 
an estimate of gravitational strength at this point. The fact that 
Mcok/Mj < 1 in the fast contracting isothermal phase is the rea- 
son that the collapsing cloud has not yet broken up into several 
fragments. 

With the approximation of the radial density profile Eq. <20t 
the total mass within a sphere of radius r ^ rcore is given by 



M(r) 



c 
G 

IM, 



3 (2^) ( 



(23) 



3.3 X 1016 cm 



where we used = 6.5 for the numerical example. Comparing the 
core mass Mcoic to total mass of the sphere, M, we get 



M 



1 
3 

0.08 



1/2 



pcoie 

Po 



pcoie 

Po 

1/2 



-1/2 



(24) 



Using Eq. <23> we can compute the mass accretion at different 



radii: 



M(r) 



c 
G 



(25) 



where Vr is the radial infall velocity. The gravitational acceleration 



outside the core region is g 



(Cc 0c) /t and the radial veloc- 



ity Vr ~ (? (^c 4>'c) (4vr G pcom)~^^^ /r assuming that the velocity 
changes on the dynamical time of the core region ^. Together with 
the size of the core region Eq. ilW we get a good estimate of the 
mass accretion onto the core: 



6.1 X 10" 



,3/2 



(26) 



3/2 



This assumption predicts a r"^ dependence of the radial velocity outside 
the core region. Our simulations show only a weak dependence of the ve- 
locity on the radius. Using the local dynamical time gives a velocity which 
is independent of r. From our simulations we conclude that the increase of 
the infall velocity occurs on a time scale somewhat in between the dynam- 
ical time and the local free fall time. Using either time scale give the same 
velocity and mass accretion at the core edge. 



Figure 8. Spherical averaged radial profiles of the mass accretion in the 
isothermal regime corresponding to Fig.|3] The maximum mass accretion 
occurs roughly at the core radius rcore which moves steadily towards the 
center (from run Al). 



Eq. <26> shows that the mass accretion of a collapsing Bonnor- 
Ebert-Sphere is independent of its initial mass and radius. This re- 
sult is applicable to non-merging, i.e. low mass , cloud cores and 
similar to the case of a collapsing SIS 

Although Eq. <26> suggests that the mass accretion is con- 
stant in time we observe a slowly increasing M which is due to the 
steady increase of the radial velocity. At a given radius r ^ rcore in 
the envelope the mass accretion approaches a constant value. The 
mass accretion has its maximum at the at the core radius and drops 
quickly towards the core center as the radial velocity sharply at 
radii r < rcore Fig.|8|shows the evolution of the mass accretion in 
the isothermal regime. 



3.3 Collapse with cooling: post-isothermal phase 



When the central density exceeds the critical density ricrit ~ 
10^'^ cm~^ thermal energy produced by gravitational contraction 
is no longer radiated efficiently and the core starts to heat up. The 
resulting increase of the thermal pressure slows down the contrac- 
tion of the central gas region and the effective equation of state 
becomes stiffen This is reflected in the increase of the effective 
equation of state, 7cff. Our simulations (cf. Fig.0 show a steep in- 
crease of 7eff from 1.0 to ~ 1.3 after the core density rises above 
the critical cooling density. Due to the complex cooling process 7eff 
does not stay at a fixed value but varies with time with a trend to 



decrease in the regime 10 



^ ricore ^ 10^^ cm~^. The strong vari- 
ability of 7 in Fig. Q reflects also the multiple shock occurrence 
during the collapse (cf. Fig. I19> . Note also that the data points in 
this plot do not correspond to the time steps in the simulation but to 
the larger time steps at which we got output data files. Sampling the 
data points on a higher time resolution might result in a smoother 
graph. 

In the post-isothermal phase, the density increase in the center 
of the gas cloud is dictated by the cooling time tcooi as opposed 
to the dynamical time tff. Fig.|9|shows the evolution time scale of 
the core region, i.e. ricore/ricore as a function of the density. After 
the core density reached ~ 10** cm^'' the collapse proceeds on a 
time scale that is determined by cooling process. At densities of 
> 10^" cm""* the core density increases on a nearly constant time 
scale much larger than the local free fall time m arking a phase of 
'hydrostatic equilibrium' (cf. also lLarsonll969t) . At later times the 
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Figure 9. Evolution time scale, ievol = "core/»i-coic, as a function of the core 

— 1/2 

density ricorc- The solid line shows the local free-fall time, i.e. t oc n^oie 
(from run Al). 

temperature of the core (at n> 10^'^ cm~'^ and 250 K) rises 
to a point when cooling by molecular excitations becomes more 
efficient (due to the abundant production of H2O and its multitude 
of excitation levels) and the collapse proceeds again almost on the 
local free fall time. 

The slower contraction of the warmer core region, as com- 
pared to the isothermal case, causes the supersonic infalling mate- 
rial to shock at the edge of the warm core. These first shocks oc- 
cur preferentially above and below the disk plane where the gas 
is not rotationally supported at a typical distance of ~ 650 AU 
(for the runs An) from the center. Temperatures at this point reach 
70 — 80 K and the temperature profile is discontinuous rising from 
30 K to 70 K between the outer pre-shock region and the shocked 
region within a few tens of AU (cf. Fig. llO> . The velocities at this 
point reach values of 1.5 km s^^ (A4 ~ 2.5). As infalling gas piles 
up at the shock boundary, a density discontinuity also develops in 
the post-shock region which separates a low density envelope and 
a high density pre-protostellar disk. This post-shock region slowly 
moves towards the center of the gas cloud and its temperature in- 
creases further while supersonic gas from the outside envelope hits 
this shock boundary. 

Gas inside this first shock region is subsonic and therefore 
not shock heated leading to a lower temperature region behind the 
shock. As the core density and the core temperature continue to in- 
crease, gas inside the pre-protostar region becomes supersonic and 
a secondary shock develops at ~ 100 AU above and below the disk 
plane (cf. second panel of Fig. fTol . Typical temperatures at the sec- 
ondary inner shock boundaries are ^ 150 K whereas the core stays 
cooler. This stage (the core density is now ^ 10^^ cm"'') marks the 
first development of the protostellar disk. 

Similar to the outer first shock, the inner shock fronts move 
slowly towards the disk plane while the local Mach number rises 
to ~ 4 corresponding to a shock velocity of 3kms"^. Again, 
the collapsing core region is fed only by gas inside the inner shock 
fronts whereas gas falling from outside onto the shocks is stalled at 
the shock boundaries and heats up the gas in these shock regions. 

Fig. \H\ shows the time evolution of the density profile and 
temperature profile, respectively, in the non-isothermal regime 
from simulation Al. The second line from the bottom in the den- 
sity plot marks the last stage of the isothermal collapse phase. In the 
non-isothermal phase the density profile deviates from the pro- 
file because of the slower contraction of the core during this regime. 
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Figure 10. Temperature line profiles along the z-axis through the center of 
the gas cloud for run Al. The first panel shows the line profile after the 
first shock develops when the central density reached 2 X 10^ cm~^. The 
second panel shows the temperature profile at a later stage when Wcoic ^ 
10^^ cm""^. At this time a second shock in the inner core region develops 
when supersonic material falls onto the protoplanetary disk (from run Al). 

The density falls of more steeply than r^'^ with increasing radius in 
the heated up region and approaches a r"^ profile in the isothermal 
regime at r > 10"^" cm. The transition region appears at a radius 
of ~ 7 X lO'^^ cm which is in good agreement with the theoretical 
prediction of Eq. <2H . The temperate profile rises steeply with de- 
creasing radius, roughly as r"^, whereas the core region has a flat 
temperature profile. The profiles p(r) are radially bin-ed spherical 
averages, i.e. p{r) = Jdf2p(x), where dfl = dcos6'd(/3 is the 
differential solid angle. 

Fig. ll2l shows the evolution of radial infall velocity as a func- 
tion of radius from run Al. The peak velocity increases with time 
and 'moves' towards smaller radii (the first line from the top refers 
to the second line from the bottom of Fig. llll which marks the end 
of the isothermal regime). Due to the shock which develops at the 
edge of the core the velocity profile becomes steeper with time 
whereas the peak velocity approaches a constant value of ~ 3 c 
in the spherical collapse phase. The velocity at the center drops al- 
ways to zero. 

In Fig. lI3l we show the evolution of the ID density profiles for 
run B68. Because the initial core density po in this case is higher 
than for the runs Al - A3 (pBes/pAi ~ 300) the isothermal region 
with p oc is smaller but the kink in the profile which sepa- 
rates the non-isothermal core region from the isothermal envelope 
appears at the same physical radius of 7 x 10^^ cm. This result is 
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Figure 11. Radial profiles of the density (upper panel) and temperature 
(lower panel) at different times. As the simulation resolution increases with 
time due to the Jeans refinement criterion the shown line profiles also in- 
crease in resolution. The dashed short line shows a density profile. The 
data are compiled from run Al. 
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Figure 12. Radial profiles of the radial velocity at different times (this pro- 
files correspond the those shown in Fig. II II . During the first stages of the 
spherical collapse the radial velocity increases steadily with time until it 
reaches the maximum value of ~ 3 c and the peak velocity is moving to- 
wards the inner part of the cloud. After the temperature in core region begins 
to increase the maximum of the radial Mach number stays at r ~ 10^*^ cm 
(from run Al. 




r [cm] 



Figure 13. Similar to figures Fig. II H and Fig. ll2l exceDt that the profiles are 
for run B68. The dashed line in the density panel shows a profile. 



again in agreement with the theoretical prediction. Finally, we note 
that a similar double shock struct ure has been seen in 2D collapse 
simulations bv lYorke et alJ il995h that feature dust as the coolant. 



4 DISK FORMATION 

As described in the previous section, the appearance of the sec- 
ondary shock establishes the initial state for the formation of the 
protostellar disk. The protostellar disk builds up and achieves a col- 
umn density of E ^ 10^'^ gcm"^ and a disk height h ~ 200 AU. 
As the shock fronts move towards the equatorial plane, the density 
and temperature increase on a time scale slower than the local free 
fall time. Fig. 114! (run Al) and Fig. I17l (run B68) show the time 
evolution of the column density which we compute as follows: 

E{R)= [ dzp{x,y,z) . (27) 



where R — y^x^ + y^ is cylindrical radius and we evaluate the 
integral <27t throughout the entire simulation box. Similar to the 
radial density profile the column density profile develops a enve- 
lope and a flat core region. The non-isothermal core is separated 
from the envelope by a steep increase in density which leads to a 
kink in the radial density profile. In general, we find that the col- 
umn density falls off roughly as r~^. The disk height decreases 
while the disk gets denser and reaches a value of ~ 10 AU when 
the column density becomes ~ 10'' g cm~^. 
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Figure 14. Radial profiles of the column density S and disk height h at 
different times. For comparison, the dashed line shows a profile. The 
shown profiles correspond to the density shown in Fig. II IK run Al). 



In the case of run B68 the ring structure is clearly visible in 
latest density profile. 

Once a disk-like object forms, one can define the disk height 

hhy 



h{R) = 



(28) 



Fig. 1141 shows the evolution of h for run Al with time. The ini- 
tial thick protostellar disk flattens to solar nebular disk with a disk 
height of a few tens of AU. 

At the time when the protostellar disk first forms, it rotates 
close to a solid body with a angular velocity of f2 ~ 10~^^ rad s~^ 
whereas the outer envelope follows a rotation law according to its 
density distribution, i.e Q oc r~^. Fig. llSl shows the evolution of 
the angular velocity for run Al: while the disk is slowly accret- 
ing gas to a core density of p ~ lO^^^gcm^^ its core spins 
up to n ~ lO^'rads^^. These numbers show that the core re- 
gion is rotationally supported against gravitational collapse, i.e. 
n ~ y/4nGp/3. 

We show the mass accretion rate of the disk, M, in Fig. ll6l (run 
Al) and Fig. ll7l (run B68). This quantity is defined by the standard 
vertically averaged continuity equation: 

dM 



(29) 



where Vr is the radial velocity in the disk. In the case of run Al 
after the protostellar disk forms the mass accretion at the outer edge 




Figure 15. Radial profiles of the angular velocity Q, the toroidal velocity 
Vy,, and the specific angular momentum at different times. The shown 
profiles correspond to the density shown in Fig. ll4l (run Al). 



of the disk is ~ 5 x 10"'' Mq /year and increases slightly with 
time to ^ 10~^ Mq /year while the disk radius decreases. For run 
B68 the mass accretion is highest at the ring's outer radius and 
reaches a maximal value of ~ 3 x 10"* AfQ/year there. Most of 
the infalling material is deposited at this outer rim of the disk and 
the mass accretion decreases towards the disk center as M oc r. 
The infalling gas outside of the disk reaches a constant terminal 
velocity of the order of the sound speed (cf. the radial velocity plot 
of Fig. I12t which results in a constant mass accretion in time in 
this outer region. In the second panel of Fig. I16l we compile the 
evolution of the total mass in the disk, A/di.sk(i?) = 27r/di?i?E, 
for run Al which reaches a few solar masses. 
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Figure 16. Radial profiles of the mass accretion M and the cumulative disk 
mass Mjisij at different times. The shown profiles correspond to the density 
shown in Fig. ll IK run Al). 

5 DISK EVOLUTION: BARS, RINGS, AND 
FRAGMENTATION 

Although we started all runs with a 10% m = 2 density perturba- 
tion, we did not observe any signs of fragmentation or non-central 
collapse in the isothermal regime. This is due to low mass (com- 
pared to the Jeans mass) in the collapsing core region (cf. Sec. l3.2> . 
The first features which appear are a bar-type or ring-type struc- 
ture. They are first evident after that stage of disk formation during 
which the non-isothermal core slowly accretes the surrounding gas. 

It is known that t hin disks are unstable to fragmentation if its 
Toomre Q -parameter jToomrell964) 



becomes smaller than unity, where k is the epicyclic frequency: 

The first unstable modes which are expected to grow and pos- 
sibly persist for some dynamical times are the m = ring-mode 
and m = 2 bar-mode. We find that a ring-type structure develops in 
the case with tff f2 = 0.2 (runs Al and B68) with a dimensionless 
radius of ^,ins ~ 6 x 10~^. For run A2 this corresponds to a physi- 
cal radius of ~ 300 AU and for run B68 to ~ 18 AU. In both cases 
the ring- structure persists only a few dynamical times and either 
fragments into a binary system (run A2) or collapses to a bar (run 
B68). In the case with tffQ, = 0.1 (run Al) and ts^l = 0.3 (run 
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Figure 17. Similar figures than Fig. ll6l and Fig. ll4l for run B68. The ring 
structure is clearly visible in the surface density at late stages. The dotted 
line shows a profile. 



A3) no such ring structures develop in the simulations. Instead a 
bar-like structure with spiral arms forms without a intermediate for- 
mation of a ring structure. A ri ng structure in the Coalsack Globule 
2 was recently discovered by iLada et alj J2004'). Our simulations 
may pertain to this system; in particular, since this cloud can be 
well described with a sub-critical Bonnor-Ebert profile. 

We show 2D images of the density, temperature, and Q- 
parameter from the runs An and B68 at different times in Figs. 1181 
- 1291 . The bar that forms in the slow rotating run Al shows no 
signs of fragmentation although the Q-parameter drops below unity 
when the bar density is very high, > 10~** g cm~^. Slight instabil- 
ities in this high density and temperature (600 — 800 K) regime 
cannot grow quickly enough to lead to fragmentation as the ineffi- 
cient cooling prevents a fast collapse of overdensities. In the case 
of faster rotation (run A2), instabilities form earlier resulting in a 
smaller Q-value than in run Al . These instabilities are large enough 
to lead to the fragmentation of the ring-type structure. In this situ- 
ation, we observe the break-up into two fragments where each of 
them has a disk with spiral structure. These disks are surrounded 
with a circumstellar disk, and are connected by a tidal stream. The 
separation between the members of this binary system is ~ 260 AU 
which is roughly twice the size of the spiral structures itself. Each 
of these protostellar systems has a total mass of ~ 1.3 A/q. Run 
A3 with tffQ, — 0.3 (Figs. l24l - I26> does not form a ring structure 
but collapses to a very twisted spiral arms and shows no sign of 
fragmentation. The simulation of Barnard 68 (run B68, Figs. 1271 - 
I29t with iff f2 — 0.2 results in a shortly lived ring structure which 
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collapses to a bar. This bar has a size of ~ 130 AU and a total mass 
of ~ 0.1 Mq. 

The double-shock structure of these disks is easily seen in 
the density plots, as well as the temperature structure, and veloc- 
ity fields in the x-z plane cuts. 

The complex temperature structure due to cooling and the re- 
sulting pressure response in the non-isothermal regime prevents a 
simple analysis to decide which kind of structure one can expect 
from a rotating collapsing cloud. From our simulations we con- 
clude that during the formation of the pre-stellar disk a ring struc- 
ture appears at a radius of 



(32) 



whenever the cloud rotates fast enough. In order to decide whether 
such possibly formed ring will fragment into several pieces, one 
has to compare this radius with the size of the non-isothermal core 
where the enhanced thermal pressure stabilizes the core region 
against fragmentation. In section ITTI we showed that the core ra- 
dius of a collapsing Bonnor-Ebert-Sphere is (almost) independent 
of the initial cloud parameters and only determined by the density 
where tff ~ tcooi (cf. Eq. I2U . For a ring to fragment its size must not 
be much smaller than this warm core size which gives the condition 



PO 
pcore 



1/2 



1/2 



(33) 



We find that ring structures which do not obey the condi- 
tion <33> collapse to a single bar rather than fragmenting into two 
or more pieces. Fies. l27l to l28l show the time evolution of the den- 
sity and temperature of the simulation B68 where a firstly formed 
ring condense to a bar. Note also that the core region rotates with 
Q{r) ~ const, whereas the disk rotates differentially outside the 
core region (cf. Fig. I15> where shearing effects enhance instabil- 
ities. Another mechanism which facilitates the fragmentation of a 
ring structure is the appearance of a shock at the core edge. A ring 
structure of similar size as the core can accrete more material than 
smaller rings increasing the possibility of fragmentation. 

In the case where tufl = 0.1 (run Al), we find that a 2- 
armed spiral structure develops when pcorc ~ lO^^^gcm^'^ and 
T ~ 200 K. If the initial sphere rotates with tffQ = 0.2 (run 
A2) a ring structure forms earlier when the density and temper- 
ature are respectively ^ lO^^^gcm"^ and ~ 85 K. This con- 
firms the general picture that slowly rotating clouds collapse to bars 
and thin filaments whereas faster rotating objects form a ring-type 
structure. Recent numerical investigation of rotating Bonnor-Ebert- 
Spheres bv lMatsumoto & Hanawal72003h showed that spheres with 
tff fl ~ 0.2 collapse to a ring which fragments into several pieces. 

The theory of b ars in disk-like structures (see 
iBinnev & Tremain3 Il987t) predicts that the bar velocity f2b 
must be larger than the so called pattern velocity: 

= n - i K (34) 

of the disk. We show the evolution of the pattern velocity as a func- 
tion of disk radius in Fig. 1301 As the bar size, Rb is limited within 
the corotation radius Rco, where the corotation radius is given by 
^l{Rco) ~ f^b, and the disk velocity falls off at the outer core the 
bar velocity is limited by: 



rip < f^b < n 



(35) 



The disk velocity in the core region is roughly independent of 
the radius r which results in a small pattern velocity, e.g. ilp <^ il. 
Due to the different time scales of the core evolution ~ tff and the 
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Table 2. Summary of the simulation results. Where the disk profiles are 
measured in the non-isothermal regime, is the bar velocity, and M* is 
the central mass of the object(s). 
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Figure 30. Evolution of the pattern velocity Qp as a function of radius R. 
The data are taken from ran Al . 



rotation period ~ Q^^ where tff S7 <C 1 we were not able to fol- 
low a full rotation of the condensed bars. Nevertheless, we estimate 
the average bar speed for our different simulations (see table|2j as, 
few X 10""' (run Al), few x lO'^^ (run A3), and ~ 10"^" rad s''^ 
(run B68) which is in agreement with condition <35> . Note that the 
run with the fastest initial rotation (run A3, tff SI = 0.3) results in 
the slowest bar velocity. This might be due to the fact the the bar 
of run A3 is much larger than the bar of run Al (tff fl — 0.1) and 
its pronounced spiral arms released sufficient angular momentum 
to the surrounding disk. 

Finally, we note that the column density profiles of the disks 
formed in our simulations are rather steep, obeying 



E oc J?~ 



(36) 



This is steeper than Hayashi models jHavashi|[l98lll but in good 
agreement with disk models inferred from the measurement by 
iKuchnei) i2004 of planets in extrasolar systems, who found E oc 

jj-2.0±0.5 



6 SUMMARY AND DISCUSSION 

In this work we studied the collapse of rotating marginal sta- 
ble Bonnor-Ebert-Spheres using 3D hydrodynamical simulations 
based on AMR technique. Compared to former numerical studies 
of collapsing molecular clouds, we included the effect of radiative 
cooling by molecular line emissions. This more realistic approach 
shows that the effective equation of state is a complex function of 
density and time during the collapsing phase and can not be approx- 
imated by a time independent EOS. Initially the gas cloud collapses 
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Figure 18. Shows 2D slices of the mass density (logarithmic scale in [gcm~^]) and velocity field at different times for run Al. After the protostellar disk 
forms a spiral feature builds up and collapses to a thin filament. The left panels and right panels show the evolution in the disk plane and perpendicular to the 
disk plane, respectively. Note that the highest resolution areas correspond (from top to bottom) to 256, 512, and 1024 pixels in each dimension. 




Figure 19. Evolution of the temperature field (gray scale in [K]) and density (contour lines). The panels correspond to these of Fig. llSK run Al). 
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Figure 20. Shows the evolution of the Toomre Q-parameter (left panel, logarithmic scale) and the angular velocity Q (right panel, logarithmic scale in 
[rads~^]. The panels correspond to these of Fig. fT8l (run Al). 
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Figure 21. Shows 2D slices of the mass density (logarithmic scale in [g cm ^]) and velocity field at different times for run A2. 
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Figure 22. Evolution of the temperature field (gray scale in [K] ) and density (contour lines). The panels con'espond to these of Fig. l21l (mn A2). 
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Figure 23. Shows the evolution of the Toomre Q-parameter (left panel, logarithmic scale) and the angular velocity Q (right panel, logarithmic scale in 
[rads~^]. The panels correspond to these of Fig. l2T1 rrun A2). 
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Figure 24. Shows 2D slices of the mass density (logarithmic scale in [g cm ^]) and velocity field at different times for run A3. 
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Figure 25. Evolution of the temperature field (gray scale in [K] ) and density (contour lines). The panels con'espond to these of Fig. l24l (mn A3). 
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Figure 26. Shows the evolution of the Toomre Q-parameter (left panel, logarithmic scale) and the angular velocity f2 (right panel, logarithmic scale in 
[rads~^]. The panels correspond to these of Fig.|24](run A3). 
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Figure 27. Shows 2D slices of the mass density (logarithmic scale in [gem '']) and velocity field from the simulation B68 at different times: from top to 
bottom: t = 5.35 , 5.40 , 5.41 Uf corresponding 3.597 , 3.628 , 3.633 x 10^ years. 
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Figure 28. Evolution of the temperature field (gray scale in [K]) and density (contour lines). The panels con'espond to these ot"Fig. l27l (run B68). 
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Figure 29. Shows the evolution of the Toomre Q-parameter (left panel, logarithmic scale) and the angular velocity Q (right panel, logarithmic scale 
[rads~^]. The panels correspond to these of Fig. lTTl rrun B68). 
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isothermally on a free fall time until the core density reaches the 
critical point where cooling becomes less efficient (i.e. tcooi > tss) at 
n ~ 10^'^ cm^''. During this isothermal phase the cloud collapses 
from outside-in where the density profile approaches to a (steadily 
increasing) flat core region with a r^^ envelope. The infall velocity 
peaks at the edge of the core where it reaches ~ 3 and drops 
to zero at the center of the gas cloud. Subsequently the tempera- 
ture in core region starts to rise and an increasing pressure support 
stabilizes the gas clump preventing immediate fragmentation. As 
cold gas from the envelope falls onto the warm core region shock 
fronts build up, separating the cold envelope from the heated up 
core. All our simulations show a similar double shock structure: an 
early outer shock and a later inner shock. Both shock fronts move 
slowly toward the core center where the inner shock region sets the 
conditions for the first appearance of the protostellar disk. The ap- 
pearance of these shocks determine also the accretion rate onto the 
core (resp. bar) as material outside the shock region is stalled at 
the shock boundaries. We find peak accretion rates for the low and 
high mass systems of 3 x 10~* MQ/year and ~ 10~^ MQ/year, 
respectively. In this work we did not keep track of the composi- 
tion of the gas during the collapse of the molecular cloud, but its is 
known tha t the chemistry might be altered within the shock region 
(e.g. Jeir gensen et ali20o4) . 

Depending on the initial rotation and density a ring mode or 
a bar develops in the disk plane. Only one of our four simulations 
result into the fragmentation of the protostellar disk. In this partic- 
ular case, a ring structure developed first which then fragments into 
a binary system in which each system contains a bar structure. The 
preferred appearance of bars in our simulations supports theoretical 
predictions that the most likely growing instability mode in disks is 
the m = 2 bar mode. 

Due to the limitations of our numerical scheme we are not able 
to follow the bar evolution for several rotation periods as the time 
scales of the bar rotation and core evolution diverge (i.e. t,ot 3> iff) 
shortly after the bar forms. Further fragmentation of the bars might 
happen during the phase of hydrogen dissociation (T > 2000 K) 
where the effective equation of state drops below 4/3 or when tidal 
interactions of the spiral amrs become stronger. 

A comparison of our A2 simuations - which produces a dis- 
tinct ring that fragments into a bin ary - with the recent observation 
of a massive disk in M 17 lofini et alJI2004D is interesting. The 
observations reveal an inner torus (within a 100 Mq disk) whose 
radius is ~ 3.8 x 10^^ cm. The radius of the dense inital ring that 
forms in our simulation A2 is of the same order (~ 10^® cm, see 
Fig. l2U which is in good agreement with these observations. 



ACKNOWLEDGEMENT 

The authors would like to thank Tom Abel and James Wadsley for 
useful and inspiring discussions. We are grateful to David Neufeld 
and his collaborators for providing us with their cooling data. 
Thanks also to the referee, Andreas Burkert, for useful comments. 
The FLASH code was in part developed by the DOE-supported 
Alliances Center for Astrophysical Thermonuclear Flashes (ASCI) 
at the University of Chicago. Our simulations were carried out on 
a 128 CPU AlphaServer SC, which is the McMaster University 
node of the SHARCNET HPC Consortium. R.B. is supported by 
the SHARCNET Postdoctoral Fellowship programm, and R.E.P. 
is supported by the Natural Sciences and Engineering Research 
Council of Canada. 



REFERENCES 

Alves J. R, Lada C. J., Lada E. A., 2001, Nature, 409, 159 

Bate M. R., Burkert A., 1997, MNRAS, 288, 1060 

Binney J., Tremaine S., 1987, Galactic dynamics. Princeton, NJ, 

Princeton University Press, 1987, 747 p. 
Bodenheimer P., Burkert A., Klein R. I., Boss A. P., 2000, Proto- 

stars and Planets IV, 675 
Bodenheimer P, Sweigart A., 1968, ApJ, 152, 515 
Bonnor W. B., 1956, MNRAS, 116, 351 
Boss A. P, 1993, ApJ, 410, 157 
Burkert A., Bodenheimer P, 1993, MNRAS, 264, 798 
Chandrasekhar S., 1967, An introduction to the study of stellar 

structure. New York: Dover, 1967 
Chini R., Hoffmeister V., Kimeswenger S., Nielbock M., 

Nurnberger D., Schmidtobreick L., Sterzik M., 2004, Nature, 

429, 155 
EbertR., 1955, ZAp, 37,217 
Foster P N., Chevalier R. A., 1993, ApJ, 416, 303 
Fryxell B., Olson K., Ricker P., Timmes F. X., Zingale M., Lamb 

D. Q., MacNeice P, Rosner R., Truran J. W., Tufo H., 2000, 
ApJS, 131,273 

Harvey D. W. A., Wilner D. J., Lada C. J., Myers PC, Alves J. R, 

Chen H., 2001, ApJ, 563, 903 
Hayashi C, 1981, Progress of Theoretical Physics Supplement, 

70, 35 

Hayashi C, Nakano T., 1965, Progr. Theoret. Phys., 34, 754 

Hunter C, 1977, ApJ, 218, 834 

Inutsuka S., Miyama S. M., 1992, ApJ, 388, 392 

J0rgensen J. K., Hogerheijde M. R., Blake G. A., van Dishoeck 

E. R, Mundy L. G., Schoier R L., 2004, A&A, 415, 1021 
Kuchner M. J., 2004, astro-ph/0405536, to appear in ApJ 

Lada C. J., Bergin E. A., Alves J. R, Huard T. L., 2003, ApJ, 586, 
286 

Lada C. J., Huard T. L., Crews L. J., Alves J. P., 2004, astro- 

ph/0404054, to appear in ApJ 
Larson R. B., 1969, MNRAS, 145, 271 
— , 2003, Reports of Progress in Physics, 66, 1651 
Matsumoto T., Hanawa T., 2003, ApJ, 595, 913 
McLaughlin D. E., Pudritz R. E., 1996, ApJ, 469, 194 
Neufeld D. A., Kaufman M. J., 1993, ApJ, 418, 263 
Neufeld D. A., Lepp S., Melnick G. J., 1995, ApJS, 100, 132 
Olson K. M., MacNeice P, Fryxell B., Ricker P, Timmes F. X., 
Zingale M., 1999, Bulletin of the American Astronomical Soci- 
ety, 31, 1430 
Penston M. V., 1969, MNRAS, 145, 457 
Racca G., Gomez M., Kenyon S. J., 2002, AJ, 124, 2178 
Shu R H., 1977, ApJ, 214, 488 

Tilley D. A., Pudritz R. E., 2004, astro-ph/0404054, to appear in 

MNRAS 
Toomre A., 1964, ApJ, 139, 1217 

Truelove J. K., Klein R. I., McKee C. R, HoUiman J. H., Howell 

L. H., Greenough J. A., 1997, ApJ, 489, L179-H 
Truelove J. K., Klein R. I., McKee C. R, HoUiman J. H., Howell 

L. H., Greenough J. A., Woods D. T, 1998, ApJ, 495, 821 
Yorke H. W., Bodenheimer P, Laughlin G., 1995, ApJ, 443, 199 
Yorke H. W., Sonnhalter C, 2002, ApJ, 569, 846 



